On the Distribution of Plasmoids In High-Lundquist-Number Magnetic Reconnection 
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The distribution function f(ip) of magnetic flux ip in plasmoids formed in high-Lundquist-number 
current sheets is studied by means of an analytic phenomenological model and direct numerical 
simulations. The distribution function is shown to follow a power law f(ip) ~ •> which differs from 
other recent theoretical predictions. Physical explanations are given for the discrepant predictions 
of other theoretical models. 



In recent years, significant advances have been made 
in understanding the role of plasmoids (or secondary is- 
lands) in magnetic reconnection, which is believed to be 
the underlying mechanism of energy release for phenom- 
ena such as solar flares, magnetospheric substorms, and 
sawtooth crashes in fusion plasmas[lj. Plasmoids often 
i form spontaneously in resistive magnetohydrodynamics 
(MHD) 0-SljHall MHD 0,0, and kinetic particle-in- 
cell (PIC) |8|4lQ( simulations of large scale reconnection. 
Evidences of plasmoids have also been found in the mag- 
netotail and the solar atmosphere [1 lMl 3] . where they 
are demonstrated to play a significant role in particle 
acceleration[3| • 

In the framework of resistive MHD, magnetic reconnec- 
tion is governed by the Lundquist number S = VaL/t], 
where Va is the upstream Alfven speed, L is the re- 
connection layer length, and rj is the resistivity. The 
classical Sweet-Parker theory (lBL flr3 | assumes the ex- 
istence of a stable, elongated current sheet and yields 
the reconnection rate ~ BVa/VS, where B is the up- 
stream magnetic field. However, it has been shown re- 
cently that when S is above a critical value S c ~ 10 4 , 
the Sweet-Parker current sheet becomes unstable to the 
plasmoid instability, with a growth rate that increases 
with 50, The reconnection layer changes to a chain 
of plasmoids connected by secondary current sheets that, 
in turn, may become unstable again. Eventually the re- 
connection layer will tend to a statistical steady state 
characterized by a hierarchical structure of plasmoids 
[l8j |. Scaling laws of the number of plasmoids n p , the 
widths S and lengths I of secondary current sheets have 
been deduced from numerical simulations. These scal- 
ing laws can be understood by noting that the process of 
break-up of the secondary current sheet will stop when 
the local Lundquist number of a secondary current sheet 
drops below S c . Assuming that all secondary current 
sheets are close to marginal stability, it can be deduced 
that I ~ tjSc/Va ~ LS C /S, 5 ~ l/y/s~ c ~ LSl /2 /S, and 
n p ~ L/l r^> 5/5 c . The reconnection rate may be es- 
timated as r] J ~ rjB/S ~ BVaI\]~S~c , independent of 

The discovery of the surprising scaling properties of 



the plasmoid instability in the linear as well as nonlin- 
ear regimes, and the ubiquity of the instability in colli- 
sional as well as collisionless regimes have raised inter- 
est in seeking a statistical description of the plasmoid 
dynamics in recent literature |20H23l|. However, exist- 
ing theoretical models give conflicting predictions. Using 
a heuristic argument based on self-similarity, Uzdensky 
et at. suggested that the distribution function f(ip) of 
plasmoids in terms of their magnetic fluxes ip follows a 
f(ip) ~ ip~ 2 power law [2l|. On the other hand, the ki- 
netic model of Fermo et al. [13] predicts a distribution 
function that decays exponentially in the tail. In this Let- 
ter, we employ both kinetic models and direct numerical 
simulations (DNS) of resistive MHD equations to study 
the distribution of plasmoids. We first recast the heuris- 
tic argument of Uzdensky et al in the form of a kinetic 
model, and show that its steady-state solutions exhibit 
both a f(ip) ~ ip~ 2 power-law regime and an exponential 
tail. This approach not only gives a formal derivation of 
the f(ip) ~ ip~ 2 power law, but also elucidates when the 
the power-law regime makes a transition to the exponen- 
tial tail. However, the results of DNS show a power law 
closer to f(ip) ~ ip -1 than to f(ip) ~ ip~ 2 . By careful 
analysis, we identify the physical causes for this devia- 
tion, and propose a modified kinetic equation that yields 
solutions consistent with the results of DNS. 

To fix ideas, we begin with a new model kinetic equa- 
tion for the plasmoid distribution function f(^p) as a 
function of the flux ip that yields the power-law solution 
obtained heuristically in [2l|. The distribution function 
f(ip) of the magnetic flux ip evolves in time due to the 
following four effects: (1) the fluxes of plasmoids increase 
due to reconnection in secondary current sheets; (2) new 
plasmoids are generated when secondary current sheets 
become unstable; plasmoids are lost by (3) coalescence 
and (4) by advection out of the reconnection layer. These 
effects can be encapsulated in the equation 
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Here N(ip) = f(ip f )dip f is the cumulative distribu- 
tion function, i.e. the number of plasmoids with fluxes 
larger than ip. In Eq. (pQ), the following assumptions 
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Figure 1. (Color online) The distribution function ([2]) for 
S= 10 6 , 10 8 , and 10 10 . 



Figure 2. (Color online) Plasmoid distribution functions from 
direct numerical simulations. 



have been made: (1) All secondary current sheets are 
close to marginal stability, therefore on average all plas- 
moids grow at a constant rate a ~ BVa/ 1 \fS~c- (2) When 
new plasmoids are created, they contain zero flux (repre- 
sented by the source term (^(VO? where 5(ip) is the Dirac 
S— function). (3) Plasmoids disappear upon encounter- 
ing larger plasmoids. This is represented by the loss term 
—/N/ta, where the characteristic time scale to encounter 
a larger plasmoid is estimated as ~ ta/N = L/NVa, as- 
suming the characteristic relative velocity between plas- 
moids is of the order of Va- The process of coalescence is 
assumed to be instantaneous. Note that when two plas- 
moids coalesce, the flux of the merged plasmoid is equal 
to the larger of the two original fluxes [2pj ]. Therefore, 
coalescence does not affect the value of / at the larger of 
the two fluxes. (4) Lastly, plasmoid loss due to advection 
is represented by the term — //ta, where the time scale 
ta is based on the outflow speed ~ Va- 

Under steady- state conditions, Eq. (pQ) admits the an- 
alytic solution 



2C/ar A 



(C - exp(-^/ar A )y 



■exp(-i/j/aT A ), (2) 



where the constant C = 1 + 2/n p , with the total number 
of plasmoids n p = J °° f(ip)dip. The source term C^(V0 
sets the boundary condition /(0) = (/a, which gives the 
relation (ta = n 2 /2-\-n p . The source term magnitude 
( may be estimated by the relation n p ~ S/S c . In the 
limit S > S c , we have ( ~ n 2 p /2r A ~ (S/S c ) 2 /2r A - The 
distribution function (|2j) has three distinct regimes when 
S > S c : (i) / — (2/olta) exp(— ^/ckta) when ip/arA > 
1; (ii) / ~ 2ar A ^~ 2 when 2/n p <C ip/ar A <C 1; (hi) / ~ 
rip/2arA when ifj/arA <C 2/n p . Therefore, the solution 
admits both an exponential tail and a power-law regime. 
It can be shown that the dominant loss mechanism in the 



former regime is advection (N <C l;adf/drf ^ — //ta), 
while it is coalescence in the latter (TV ^> 1; ad f /dip c± 
—JN/ta)- Figure [1] shows the distribution function (|2j) 
for 5 = 10 6 , 10 8 , and 10 10 . Here, to fix ideas, we have 
taken S c = 10 4 , Va = 1, B = 1, and L = 1. and all 
scaling relations, such as a ~ BVa/VS^, are replaced by 
equalities. Note that the range where the f ~ ip~ 2 power 
law holds is more extended for higher S. 

To test the f(ip) ~ i(j~ 2 power law by DNS, we use the 
same simulation setup of two coalescing magnetic islands 
as in a previous study [19] . The 2D simulation box is the 
domain (x,z) G [-1/2,1/2] x [-1/2,1/2]. In normalized 
units, the initial magnetic field is given by Bo = V^o x y, 
where ipo = tanh (z/h) cos (ttx) sin (2ttz) /2tt. The pa- 
rameter h, which is set to 0.01 for all simulations, deter- 
mines the initial current layer width. The initial plasma 
density p is approximately 1, and the plasma tempera- 
ture T is 3. The density profile has a weak nonuniformity 
such that the initial condition is approximately force- 
balanced. The initial peak magnetic field and Alfven 
speed are both approximately unity. The plasma beta 
f3 = p/B 2 = 2pT/B 2 is greater than 6 everywhere. Per- 
fectly conducting and free slipping boundary conditions 
are imposed along both x and z directions. Only the 
upper half of the domain {z > 0) is simulated, and so- 
lutions in the lower half are inferred by symmetries. We 
use a uniform grid along the x direction and a nonuni- 
form grid along the z direction that packs high resolu- 
tion around z = 0. For cases with S = 10 6 and 3 x 10 6 , 
the mesh size is 12726 x 1600, and the smallest grid size 
along z is 5.7 x 10 -6 . For the S = 10 7 case, the mesh 
size is 37800 x 2880, and the smallest grid size along z is 
1.9 x 10 -6 . No explicit viscosity is employed in these sim- 
ulations. A fourth order numerical dissipation is added 
to damp small fluctuations at grid scale [24] . 

The initial velocity is seeded with a random noise of 
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Figure 3. (Color online) The plasmoid distribution with re- 
spect to the relative speed Av and the flux ip from the run 
S= 10 7 . 



amplitude 10 -6 to trigger the plasmoid instability. The 
early period when the reconnected flux is less than 0.01 
is precluded from the analysis to allow the reconnec- 
tion layer to reach a statistical steady state. We take 
data during the period when the reconnected flux is be- 
tween 0.01 to 0.05, corresponding to 25% of the initial 
flux in each of the merging islands. This period roughly 
spans 6t^, insensitive to S. Snapshots are taken at inter- 
vals of O.Olr^. We identify plasmoids within the range 
x G [—0.25, 0.25] with a computer program for each snap- 
shot, which provide the dataset for further statistical 
analysis. Figure [2] shows the probability distribution 
functions for S = 10 6 , 3 x 10 6 [two runs, labeled 

as (a) and (b)], and S = 10 7 . Distribution functions are 
normalized such that J °° f(ip)dip is equal to the average 
number of plasmoids in each time slice. These numerical 
results appear to be robust and reproducible, as exem- 
plified by the two S = 3 x 10 6 runs that yield nearly 
identical distribution functions. Qualitative similarities 
between Fig. [1] and Fig. [2j especially the existence of 
three distinct regimes, are evident. However, the dis- 
tribution function in the power-law regime is closer to 
f(ijS) rsj ijj~ l instead of f(i/j) ~ i\)~ 2 . 

To understand the discrepancy between the numerical 
results and the power-law prediction of Eq. (2), we need 
to critically examine the basic assumptions that give rise 
to the f(ip) ~ power law. In the /(VO ^ regime, 
the dominant balance in Eq. (pQ) is between the plas- 
moid growth term and the loss term due to coalescence, 
i.e. adf '/d^ ~ —/N/ta- A key assumption underly- 
ing the loss term —/N/ta is that the relative speeds of 
a plasmoid with respect to neighboring plasmoids larger 
than itself are of the order of Va and are uncorrelated 
to the flux of the plasmoid. To examine this assumption 
with numerical data, we measure the relative velocity Av 



of each plasmoid at any given time with respect to the 
first larger plasmoid it will encounter by extrapolating 
the trajectories of the plasmoids with their velocities at 
that time. Note that Av is undefined for the largest plas- 
moid, or when all larger plasmoids are moving away from 
a given plasmoid. The plasmoids with Av undefined are 
disregarded in the analyses. Figure [3] shows the distri- 
bution g(ip, Av) of plasmoids with respect to i/j and Av 
from the run S = 10 7 . Here we normalize g(tp, Av) such 
that g(i/j, Av)d(Av) = 1 for better visualization. We 
can clearly see that the distribution is not uniform across 
different values of ip. The distribution covers a broader 
range of Av at smaller -0, and it becomes more concen- 
trated around Av = at larger Similar results are 
also observed in other runs. Therefore, it appears that 
the reconnection layer organizes itself spontaneously into 
a state such that large plasmoids tend to avoid coalescing 
with each other. 

How do we interpret this phenomenon? As discussed 
earlier, the flux of a plasmoid is approximately propor- 
tional to its age because all plasmoids grow approxi- 
mately at the same rate a. Consequently, a plasmoid 
can become large only if it has not encountered plas- 
moids larger than itself for an extended period of time. 
Presumably, plasmoids moving rapidly relative to their 
neighbors will encounter larger plasmoids and disappear 
easily, whereas those with small relative speeds are more 
likely to survive for a long time and become large. This 
observation motivates us to consider a distribution func- 
tion Fty, v), where v can be interpreted as the plasmoid 
velocity relative to the mean flow (which has a profile 
along the outflow direction) . The governing equation for 
F(ip, v) is written as 

dtF + a?-=(5^)h(v) , (3) 

dip T A T A 

where the function H is defined as 

r°° i r°° \v-v'\ 
H(^v)= d^j dv — F W> v ')i (4) 

JiP J-oo V A 

and h(v) is an arbitrary distribution function in velocity 
space when new plasmoids are generated. The distribu- 
tion function f(ip) can be obtained by integrating F(ip,v) 
over the velocity space. Eq. §3§ differs from Eq. (pQ) in 
the plasmoid loss term due to coalescence, where the rel- 
ative speed \v — v'\ between two plasmoids is taken into 
account in the integral operator of Eq. (j4]) . If we replace 
\v - v'\ in Eq. (@|) by Va, then Eq. (j3]) reduces to Eq. 
(PQ). Steady-state solutions of Eq. ((3]) can be obtained 
numerically. To fix ideas, we assume a Gaussian profile 
h(v) = (I/v^Va) exp(— v 2 /Va) for the arbitrary source 
function. Fig. |4] shows the resulting f(ip) for (ta = 10 6 , 
10 7 , and 10 8 . Assuming n p c± S/S c and S c ~ 10 4 , these 
solutions approximately correspond to S = 3 x 10 7 , 10 8 , 
and 3 x 10 8 , respectively. These solutions also show three 
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Figure 4. (Color online) Distribution functions from numeri- 
cal solutions of Eq. ((3} ■ 



distinct regimes as the solutions in Fig. [TJ However, the 
distribution in the intermediate power-law regime is close 
to f(ip) ~ consistent with DNS. We have tried other 
smooth h(v) profiles, and the results do not appear to be 
sensitive to the specific form of h(v), as long as h(v) cov- 
ers a broad range of v (typically of the order of Va)- 

A previous DNS study of plasmoid distribution has 
been recently carried out by Loureiro et al. (23j . where 
they claimed confirmation of the f(ip) ~ ip~ 2 distribu- 
tion. It should be pointed out that Loureiro et al. com- 
pared the f(ip) ~ ?p~ 2 prediction with simulation data in 
the large-?/? regime. If we focus on the smaller-?/? regime of 
their numerical data, the distribution appears more con- 
sistent with our finding f(i/i) ~ ?/? _1 . This flattening of 
distribution function in the smaller-?/? regime was noted 
by Loureiro et al., but no attempt was made to fit the 
smaller-?/; regime to a power law. An important question 
is: do we expect to see a power law in the large-?/; regime 
or the smaller-?/; regime? Our analytic theory reveals 
that the transition from a power-law distribution to an 
exponential tail is due to a change in the dominant loss 
mechanism from coalescence to advection, which occurs 
approximately when TV ~ 0(1). In our simulation data, 
the cumulative distribution function N(ip) drops below 
unity at ip ~ 10 ~ 3 , which is also approximately where 
the distribution function deviates from f(ijj) ~ ?/; _1 to 
a more rapid, presumably exponential, falloff. There- 
fore, this rapidly falling tail is not where a power law 
should arise. However, simulation data in the large-?/; 
regime is sufficiently uncertain that it may be difficult to 
make a clear distinction between a ip~ 2 and an exponen- 
tial falloff. Note that the exponential falloff at large ip 
is consistent with the prediction of the kinetic model of 
Fermo et al. [20] and a subsequent analysis of the flux 
transfer events (FTEs) in the magnetopause from Cluster 
[22j . Fermo et al. did not explicitly address the distribu- 



tion of smaller plasmoids. Because the coalescence term 
in their model is based on very different considerations 
and assumptions, it is not clear whether the distribution 
of smaller plasmoids will follow a power law. 

Although Eq. ([3j) is a significant improvement on Eq. 
(PQ), it does not include some important physical effects. 
Most notably, coalescence between islands is assumed to 
occur instantaneously, whereas in reality larger plasmoids 
take longer to merge, and there can be bouncing (or slosh- 
ing) between them [25|, [26| . These effects may also con- 
tribute to the distribution shown in Fig. [3l Furthermore, 
the velocity v relative to the mean flow is assumed to 
remain constant throughout the lifetime of a plasmoid, 
whereas in reality some variation is expected due to the 
complex dynamics between plasmoids. Finally, in high-^ 
regime the current sheet between two coalescing plas- 
moids can also be the source of more plasmoids. [27] . 

It should be borne in mind that our considerations are 
valid for collisional plasmas obeying the resistive MHD 
equations. In weakly collisional systems the plasmoid in- 
stability inevitably drives reconnection towards the col- 
lisionless regime p, 0, The question of the plasmoid 
distribution in the collisionless regime remains largely 
open. However, some of the key ideas in this work, such 
as the tendency of large plasmoids to avoid coalescence, 
may still be relevant. The present study is limited to 
highly idealized 2D problems where more concrete con- 
clusions can be drawn. In 3D geometry oblique tearing 
modes have been shown to play an important role 28, 29|, 
and a statistical description of such systems remains a 
great challenge. 
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